# -*- Mode: shell-script -*- 

GlobalPackage(gap.int);

#############################################################################
##
#A  integer.g                   GAP library                  Martin Schoenert
#A                                                           & Alice Niemeyer
#A                                                            & Werner Nickel
#A                                                              & Alex Wegner
##
##
#Y  Copyright (C) 2018-2021, Carnegie Mellon University
#Y  All rights reserved.  See LICENSE for details.
#Y  
#Y  This work is based on GAP version 3, with some files from version 4.  GAP is
#Y  Copyright (C) (1987--2021) by the GAP Group (www.gap-system.org).
##
##  This file contains  those  functions  that  mainly  deal  with  integers.
##
##


#############################################################################
##
#V  Integers  . . . . . . . . . . . . . . . . . . . . .  ring of the integers
#V  IntegersOps . . . . . . . . . . . . . . . . operation record for integers
##
IntegersOps := Copy( RingOps );

Integers := rec(
    isDomain                    := true,
    isRing                      := true,

    generators                  := [ 1 ],
    zero                        := 0,
    one                         := 1,
    name                        := "Integers",

    size                        := "infinity",
    isFinite                    := false,
    isCommutativeRing           := true,
    isIntegralRing              := true,
    isUniqueFactorizationRing   := true,
    isEuclideanRing             := true,
    units                       := [ -1, 1 ],

    operations                  := IntegersOps
);


#############################################################################
##
#F  IntegersOps.Ring(<elms>)  . . . . . . . . ring generated by some integers
##
IntegersOps.Ring := function ( elms )
    return Integers;
end;


#############################################################################
##
#F  IntegersOps.DefaultRing(<elms>) . . . . . . default ring of some integers
##
IntegersOps.DefaultRing := function ( elms )
    return Integers;
end;


#############################################################################
##
#F  IntegersOps.Field(<elms>) . . . . . . .  field generated by some integers
##
IntegersOps.Field := function ( elms )
    return Rationals;
end;


#############################################################################
##
#F  IntegersOps.DefaultField(<elms>) . . . . . default field of some integers
##
IntegersOps.DefaultField := function ( elms )
    return Rationals;
end;


#############################################################################
##
#F  IntegersOps.\in(<n>,<Integers>)  . . . . . . mebership test for integers
##
IntegersOps.\in := function ( n, Integers )
    return IsInt( n );
end;


#############################################################################
##
#F  IntegersOps.Random(<Integers>)  . . . . . . . . . . . . .  random integer
##
NrBitsInt := function ( n )
    local   nr, nr64;
    nr64:=[0,1,1,2,1,2,2,3,1,2,2,3,2,3,3,4,1,2,2,3,2,3,3,4,2,3,3,4,3,4,4,5,
           1,2,2,3,2,3,3,4,2,3,3,4,3,4,4,5,2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6];
    nr := 0;
    while 0 < n  do
        nr := nr + nr64[ n mod 64 + 1 ];
        n := QuoInt( n, 64 );
    od;
    return nr;
end;

IntegersOps.Random := function ( Integers )
    return NrBitsInt( RandomList( [0..2^20-1] ) ) - 10;
end;


#############################################################################
##
#F  IntegersOps.Quotient(<Integers>,<n>,<m>)  . . .  quotient of two integers
##
IntegersOps.Quotient := function ( Integers, n, m )
    local   q;
    q := QuoInt( n, m );
    if n <> q * m  then
        q := false;
    fi;
    return q;
end;


#############################################################################
##
#F  IntegersOps.StandardAssociate(<Integers>,<n>) . . . . . .  absolute value
##
IntegersOps.StandardAssociate := function ( Integers, n )
    if n < 0  then
        return -n;
    else
        return n;
    fi;
end;


#############################################################################
##
#F  IntegersOps.IsPrime(<Integers>,<n>) .  test whether an integer is a prime
##
IntegersOps.IsPrime := function ( Integers, n )
    return IsPrimeInt( n );
end;


#############################################################################
##
#F  IntegersOps.IsIrreducible(<Integers>,<n>)
##
IntegersOps.IsIrreducible := IntegersOps.IsPrime;

#############################################################################
##
#F  IntegersOps.MinimalPolynomial( <Integers>, <n> )
##
IntegersOps.MinimalPolynomial := function( Integers, n )
    return Polynomial( Rationals, MinPol( Rationals, n ) );
    end;

#############################################################################
##
#V  Primes[]  . . . . . . . . . . . . . . . . . . . . . . . .  list of primes
##
##  'Primes' is a set, i.e., sorted list, of the 168 primes less than 1000.
##
##  This is used in 'IsPrimeInt' and 'FactorsInt' to cast  out  small  primes
##  quickly.
##
Primes := [   2,  3,  5,  7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53,
     59, 61, 67, 71, 73, 79, 83, 89, 97,101,103,107,109,113,127,131,137,139,
    149,151,157,163,167,173,179,181,191,193,197,199,211,223,227,229,233,239,
    241,251,257,263,269,271,277,281,283,293,307,311,313,317,331,337,347,349,
    353,359,367,373,379,383,389,397,401,409,419,421,431,433,439,443,449,457,
    461,463,467,479,487,491,499,503,509,521,523,541,547,557,563,569,571,577,
    587,593,599,601,607,613,617,619,631,641,643,647,653,659,661,673,677,683,
    691,701,709,719,727,733,739,743,751,757,761,769,773,787,797,809,811,821,
    823,827,829,839,853,857,859,863,877,881,883,887,907,911,919,929,937,941,
    947,953,967,971,977,983,991,997 ];
IsSet( Primes );


#############################################################################
##
#V  Primes2[] . . . . . . . . . . . . . . . . . . . . . additional prime list
##
##  'Primes2' contains those primes found by  'IsPrimeInt'  that  are  not in
##  'Primes'.  'Primes2' is kept sorted, but may contain holes.
##
##  'IsPrimeInt' and  'FactorsInt'  use this list to   cast out already found
##  primes quickly.  If 'IsPrimeInt' is called only  for random integers this
##  list  would  be quite useless.   However, users  do not  behave randomly.
##  Instead, it is not uncommon to factor the  same integer twice.  Likewise,
##  once we have  tested  that $2^{31}-1$  is prime,  factoring $2^{62}-1$ is
##  very cheap, because the former divides the latter.
##
##  This list is initialized to contain all the prime factors of the integers
##  $2^n-1$ with $n < 201$,  $3^n-1$ with $n < 101$,  $5^n-1$ with $n < 101$,
##  $7^n-1$ with $n < 91$, $11^n-1$ with $n < 79$, and $13^n-1$ with $n < 37$
##  that are larger than $10^7$.
##
Primes2 := [
10047871, 10567201, 10746341, 12112549, 12128131, 12207031, 12323587,
12553493, 12865927, 13097927, 13264529, 13473433, 13821503, 13960201,
14092193, 14597959, 15216601, 15790321, 16018507, 18837001, 20381027,
20394401, 20515111, 20515909, 21207101, 21523361, 22253377, 22366891,
22996651, 23850061, 25781083, 26295457, 28325071, 28878847, 29010221,
29247661, 29423041, 29866451, 32234893, 32508061, 36855109, 41540861,
42521761, 43249589, 44975113, 47392381, 47763361, 48544121, 48912491,
49105547, 49892851, 51457561, 55527473, 56409643, 56737873, 59302051,
59361349, 59583967, 60816001, 62020897, 65628751, 69566521, 75068993,
76066181, 85280581, 93507247, 96656723, 97685839,
106431697, 107367629, 109688713, 110211473, 112901153, 119782433, 127540261,
134818753, 134927809, 136151713, 147300841, 160465489, 164511353, 177237331,
183794551, 184481113, 190295821, 190771747, 193707721, 195019441, 202029703,
206244761, 212601841, 212885833, 228511817, 231769777, 234750601, 272010961,
283763713, 297315901, 305175781, 308761441, 319020217, 359390389, 407865361,
420778751, 424256201, 432853009, 457315063, 466344409, 510810301, 515717329,
527093491, 529510939, 536903681, 540701761, 550413361, 603926681, 616318177,
632133361, 715827883, 724487149, 745988807, 815702161, 834019001, 852133201,
857643277, 879399649,
1001523179, 1036745531, 1065264019, 1106131489, 1169382127, 1390636259,
1503418321, 1527007411, 1636258751, 1644512641, 1743831169, 1824179209,
1824726041, 1826934301, 1866013003, 1990415149, 2127431041, 2147483647,
2238236249, 2316281689, 2413941289, 2481791513, 2550183799, 2576743207,
2664097031, 2767631689, 2903110321, 2931542417, 3158528101, 3173389601,
3357897971, 4011586307, 4058036683, 4278255361, 4375578271, 4562284561,
4649919401, 4698932281, 4795973261, 4885168129, 5960555749, 6809710909,
7068569257, 7151459701, 7484047069, 7685542369, 7830118297, 7866608083,
8209475377, 8831418697, 9598959833,
10879733611, 11898664849, 12447002677, 13455809771, 13564461457, 13841169553,
13971969971, 14425532687, 15085812853, 15768033143, 15888756269, 16148168401,
17154094481, 17189128703, 19707683773, 22434744889, 23140471537, 23535794707,
24127552321, 25480398173, 25829691707, 25994736109, 27669118297, 27989941729,
28086211607, 30327152671, 32952799801, 33057806959, 35532364099, 39940132241,
43872038849, 45076044553, 47072139617, 50150933101, 54410972897, 56625998353,
60726444167, 61070817601, 62983048367, 70845409351, 76831835389, 77158673929,
77192844961, 78009515593, 83960385389, 86950696619, 88959882481, 99810171997,
115868130379, 125096112091, 127522693159, 128011456717, 128653413121,
131105292137, 152587500001, 158822951431, 159248456569, 164504919713,
165768537521, 168749965921, 229890275929, 241931001601, 269089806001,
282429005041, 332207361361, 374857981681, 386478495679, 392038110671,
402011881627, 441019876741, 447600088289, 487824887233, 531968664833,
555915824341, 593554036769, 598761682261, 641625222857, 654652168021,
761838257287, 810221830361, 840139875599, 918585913061,
1030330938209, 1047623475541, 1113491139767, 1133836730401, 1273880539247,
1534179947851, 1628744948329, 1654058017289, 1759217765581, 1856458657451,
2098303812601, 2454335007529, 2481357870461, 2549755542947, 2663568851051,
2879347902817, 2932031007403, 3138426605161, 3203431780337, 3421169496361,
3740221981231, 4363953127297, 4432676798593, 4446437759531, 4534166740403,
4981857697937, 5625767248687, 6090817323763, 6493405343627, 6713103182899,
6740339310641, 7432339208719, 8090594434231, 8157179360521, 8737481256739,
8868050880709, 9361973132609, 9468940004449, 9857737155463,
10052678938039, 10979607179423, 13952598148481, 15798461357509,
18158209813151, 22125996444329, 22542470482159, 22735632934561,
23161037562937, 23792163643711, 24517014940753, 24587411156281,
28059810762433, 29078814248401, 31280679788951, 31479823396757,
33232924804801, 42272797713043, 44479210368001, 45920153384867,
49971617830801, 57583418699431, 62911130477521, 67280421310721,
70601370627701, 71316922984999, 83181652304609, 89620825374601,
110133112994711, 140737471578113, 145295143558111, 150224123975857,
204064664440913, 205367807127911, 242099935645987, 270547105429567,
303567967057423, 332584516519201, 434502978835771, 475384700124973,
520518327319589, 560088668384411, 608459012088799, 637265428480297,
643170158708221, 707179356161321, 926510094425921, 990643452963163,
1034150930241911, 1066818132868207, 1120648576818041, 1357105535093947,
1416258521793067, 1587855697992791, 1611479891519807, 1628413557556843,
1958423494433591, 2134387368610417, 2646507710984041, 2649263870814793,
2752135920929651, 2864226125209369, 4889988840047743, 5420506947192709,
6957533874046531, 9460375336977361, 9472026608675509,
12557612956332313, 13722816749522711, 14436295738510501, 18584774046020617,
18624275418445601, 20986207825565581, 21180247636732981, 22666879066355177,
27145365052629449, 46329453543600481, 50544702849929377, 59509429687890001,
60081451169922001, 70084436712553223, 76394148218203559, 77001139434480073,
79787519018560501, 96076791871613611,
133088039373662309, 144542918285300809, 145171177264407947,
153560376376050799, 166003607842448777, 177722253954175633,
196915704073465747, 316825425410373433, 341117531003194129,
380808546861411923, 489769993189671059, 538953023961943033,
581283643249112959, 617886851384381281, 625552508473588471,
645654335737185721, 646675035253258729, 658812288653553079,
768614336404564651, 862970652262943171, 909456847814334401,
1100876018364883721, 1195857367853217109, 1245576402371959291,
1795918038741070627, 2192537062271178641, 2305843009213693951,
2312581841562813841, 2461243576713869557, 2615418118891695851,
2691614274040036601, 3011347479614249131, 3358335487319458201,
3421093417510114543, 3602372010909260861, 3747607031112307667,
3999088279399464409, 4710883168879506001, 5079304643216687969,
5559917315850179173, 5782172113400990737, 6106505825833677713,
6115909044841454629, 9213624084535989031, 9520972806333758431,
10527743181888260981, 14808607715315782481, 18446744069414584321,
26831423036065352611, 32032215596496435569, 34563155350221618511,
36230454570129675721, 58523123221688392679, 82064241848634269407,
86656268566282183151, 87274497124602996457,
157571957584602258799, 162715052426691233701, 172827552198815888791,
195489390796456327201, 240031591394168814433, 344120456368919234899,
358475907408445923469, 846041103974872866961,
2519545342349331183143, 3658524738455131951223, 3976656429941438590393,
5439042183600204290159, 8198241112969626815581,
11600321878916922053491, 12812432238302009985937, 17551032119981679046729,
18489605314740987765913, 27665283091695977275201, 42437717969530394595211,
57912614113275649087721, 61654440233248340616559, 63681511996418550459487,
105293313660391861035901, 155285743288572277679887, 201487636602438195784363,
231669654363683130095909, 235169662395069356312233, 402488219476647465854701,
535347624791488552837151, 604088623657497125653141, 870035986098720987332873,
950996059627210897943351,
1412900479108654932024439, 1431185706701868962383741,
2047572230657338751575051, 2048568835297380486760231,
2741672362528725535068727, 3042645634792541312037847,
3745603812007166116831643, 4362139336229068656094783,
4805345109492315767981401, 5042939439565996049162197,
7289088383388253664437433, 8235109336690846723986161,
9680647790568589086355559, 9768997162071483134919121,
9842332430037465033595921,
11053036065049294753459639, 11735415506748076408140121,
13842607235828485645766393, 17499733663152976533452519,
26273701844015319144827917, 75582488424179347083438319,
88040095945103834627376781,
100641220283951395639601683, 140194179307171898833699259,
207617485544258392970753527, 291280009243618888211558641,
303309617049998388989376043, 354639323684545612988577649,
618970019642690137449562111, 913242407367610843676812931,
7222605228105536202757606969, 7248808599285760001152755641,
8170509011431363408568150369, 8206973609150536446402438593,
9080418348371887359375390001,
14732265321145317331353282383, 15403468930064931175264655869,
15572244900182528777225808449, 18806327041824690595747113889,
21283620033217629539178799361, 37201708625305146303973352041,
42534656091583268045915654719, 48845962828028421155731228333,
123876132205208335762278423601, 134304196845099262572814573351,
172974812463239310024750410929, 217648180992721729506406538251,
227376585863531112677002031251,
1786393878363164227858270210279, 2598696228942460402343442913969,
2643999917660728787808396988849, 3340762283952395329506327023033,
5465713352000770660547109750601,
28870194250662203210437116612769, 70722308812401674174993533367023,
78958087694609321439660131899631, 88262612316754526107621113329689,
162259276829213363391578010288127, 163537220852725398851434325720959,
177635683940025046467781066894531,
2679895157783862814690027494144991, 3754733257489862401973357979128773,
5283012903770196631383821046101707, 5457586804596062091175455674392801,
10052011757370829033540932021825161, 11419697846380955982026777206637491,
38904276017035188056372051839841219,
1914662449813727660680530326064591907, 7923871097285295625344647665764672671,
9519524151770349914726200576714027279,
10350794431055162386718619237468234569,
170141183460469231731687303715884105727,
1056836588644853738704557482552056406147,
6918082374901313855125397665325977135579,
235335702141939072378977155172505285655211,
360426336941693434048414944508078750920763,
1032670816743843860998850056278950666491537,
1461808298382111034194027645506019619578037,
79638304766856507377778616296087448490695649,
169002145064468556765676975247413756542145739,
8166146875847876762859119015147004762656450569,
18607929421228039083223253529869111644362732899,
33083146850190391025301565142735000331370209599,
138497973518827432485604572537024087153816681041,
673267426712748387612994804392183645147042355211,
1489459109360039866456940197095433721664951999121,
4884164093883941177660049098586324302977543600799,
466345922275629775763320748688970211803553256223529,
26828803997912886929710867041891989490486893845712448833,
153159805660301568024613754993807288151489686913246436306439,
1051153199500053598403188407217590190707671147285551702341089650185945215953
];
IsSet( Primes2 );


#############################################################################
##
#F  IsPrimeInt( <n> ) . . . . . . . . . . . . . . . . . . .  test for a prime
##
##  'IsPrimeInt' returns 'false'  if it can  prove that <n>  is composite and
##  'true' otherwise.  By  convention 'IsPrimeInt(0) = IsPrimeInt(1) = false'
##  and we define 'IsPrimeInt( -<n> ) = IsPrimeInt( <n> )'.
##
##  'IsPrimeInt' does trial divisions by the primes less  than 1000 to detect
##  composites with a factor less than 1000 and  primes  less  than  1000000.
##
##  'IsPrimeInt' then checks that $n$ is a strong pseudoprime to the  base 2.
##  This uses Fermats theorem which says $2^{n-1}=1$ mod $n$ for a prime $n$.
##  If $2^{n-1}\<>1$ mod $n$, $n$ is composite, 'IsPrimeInt' returns 'false'.
##  There are composite numbers for which $2^{n-1}=1$,  but they are  seldom.
##
##  Then 'IsPrimeInt' checks that $n$ is a Lucas pseudoprime for $p$, choosen
##  so that the discriminant $d=p^2/4-1$ is an  quadratic nonresidue mod $n$.
##  I.e., 'IsPrimeInt' takes the root $a = p/2+\sqrt{d}$ of $x^2 - px + 1$ in
##  the  ring $Z_n[\sqrt{d}] and computes the  traces of $a^n$ and $a^{n+1}$.
##  If $n$ is a prime, this  ring is the field of  order $n^2$ and raising to
##  the $n$th power is conjugation, so $trace(a^n)=p$ and $trace(a^{n+1})=2$.
##  However, these identities hold only for extremly few composite numbers.
##
##  Note that  this  test  for $trace(a^n) = p$  and  $trace(a^{n+1}) = 2$ is
##  usually formulated using the Lucas sequences  $U_k = (a^k-b^k)/(a-b)$ and
##  $V_k = (a^k+b^k)=trace(a^k)$, where one tests $U_{n+1} = 0, V_{n+1} = 2$.
##  However, the trace test is equivalent and requires fewer multiplications.
##  Thanks to Daniel R. Grayson (dan@symcom.math.uiuc.edu)  for  telling  me.
##
##  'IsPrimeInt' can be shown to return the correct answer for $n < 10^{13}$,
##  by testing against R.G.E. Pinch's list of all pseudoprimes to base 2 less
##  than $10^{13}$ ('ftp://cmms.cam.ac.uk/pub/rgep/PSP/psp13').
##
##  Better descriptions of the algorithm and related topics can be found  in:
##  G. Miller, cf. Algorithms and Complexity ed. Traub, AcademPr, 1976, 35-36
##  C. Pomerance et.al., Pseudoprimes to 25*10^9, MathComp 35 1980, 1003-1026
##  D. Knuth, Seminumerical Algorithms  (TACP II),  AddiWesl,  1973,  378-380
##  G. Gonnet, Heuristic Primality Testing, Maple Newsletter 4,  1989,  36-38
##  R. Baillie, S. Wagstaff, Lucas Pseudoprimes, MathComp 35 1980,  1391-1417
##  R. Pinch, Some Primality Testing Algorithms, Notic. AMS 9 1993, 1203-1210
##
TraceModQF := function ( p, k, n )
    local  trc;
    if k = 1  then
        trc := [ p, 2 ];
    elif k mod 2 = 0  then
        trc := TraceModQF( p, k/2, n );
        trc := [ (trc[1]^2 - 2) mod n, (trc[1]*trc[2] - p) mod n ];
    else
        trc := TraceModQF( p, (k+1)/2, n );
        trc := [ (trc[1]*trc[2] - p) mod n, (trc[2]^2 - 2) mod n ];
    fi;
    return trc;
end;

IsPrimeInt := function ( n )
    local  p, e, o, x, i, d;

    # make $n$ positive and handle trivial cases
    if n < 0         then n := -n;       fi;
    if n in Primes   then return true;   fi;
    if n in Primes2  then return true;   fi;
    if n <= 1000     then return false;  fi;

    # do trial divisions by the primes less than 1000
    # faster than anything fancier because $n$ mod <small int> is very fast
    for p  in Primes  do
        if n mod p = 0  then return false;  fi;
        if n < (p+1)^2  then AddSet(Primes2,n);  return true;   fi;
    od;

    # do trial division by the other known primes
    for p  in Primes2  do
        if n mod p = 0  then return false;  fi;
    od;

    # find $e$ and $o$ odd such that $n-1 = 2^e * o$
    e := 0;  o := n-1;   while o mod 2 = 0  do e := e+1;  o := o/2;  od;

    # look at the seq $2^o, 2^{2 o}, 2^{4 o}, .., 2^{2^e o}=2^{n-1}$
    x := PowerModInt( 2, o, n );
    i := 0;
    while i < e  and x <> 1  and x <> n-1  do
        x := x * x mod n;
        i := i + 1;
    od;

    # if it is not of the form $.., -1, 1, 1, ..$ then $n$ is composite
    if not (x = n-1 or (i = 0 and x = 1))  then
        return false;
    fi;

    # there are no strong pseudo-primes to base 2 smaller than 2047
    if n < 2047  then
        AddSet( Primes2, n );
        return true;
    fi;

    # make sure that $n$ is not a perfect power (especially not a square)
    if SmallestRootInt(n) < n  then
        return false;
    fi;

    # find a quadratic nonresidue $d = p^2/4-1$ mod $n$
    p := 2;  while Jacobi( p^2-4, n ) <> -1  do p := p+1;  od;

    # for a prime $n$ the trace of $(p/2+\sqrt{d})^n$ must be $p$
    # and the trace of $(p/2+\sqrt{d})^{n+1}$ must be 2
    if TraceModQF( p, n+1, n ) = [ 2, p ]  then
        AddSet( Primes2, n );
        return true;
    fi;

    # $n$ is not a prime
    return false;
end;


#############################################################################
##
#F  IsPrimePowerInt( <n> )  . . . . . . . . . . . test for a power of a prime
##
##  'IsPrimePowerInt' returns 'true' if the integer <n>  is a prime power and
##  'false' otherwise.
##
IsPrimePowerInt := function ( n )
    return IsPrimeInt( SmallestRootInt( n ) );
end;


#############################################################################
##
#F  NextPrimeInt( <n> ) . . . . . . . . . . . . . . . . . . next larger prime
##
##  'NextPrimeInt' returns the smallest prime  which is strictly larger  than
##  the integer <n>.
##
NextPrimeInt := function ( n )
    if   -3 = n             then n := -2;
    elif -3 < n  and n < 2  then n :=  2;
    elif n mod 2 = 0        then n := n+1;
    else                         n := n+2;
    fi;
    while not IsPrimeInt(n)  do
        if n mod 6 = 1  then n := n+4;
        else                 n := n+2;
        fi;
    od;
    return n;
end;


#############################################################################
##
#F  PrevPrimeInt( <n> ) . . . . . . . . . . . . . . .  previous smaller prime
##
##  'PrevPrimeInt' returns the largest prime  which is strictly smaller  than
##  the integer <n>.
##
PrevPrimeInt := function ( n )
    if    3 = n             then n :=  2;
    elif -2 < n  and n < 3  then n := -2;
    elif n mod 2 = 0        then n := n-1;
    else                         n := n-2;
    fi;
    while not IsPrimeInt(n)  do
        if n mod 6 = 5  then n := n-4;
        else                 n := n-2;
        fi;
    od;
    return n;
end;


#############################################################################
##
#F  IntegersOps.Factors(<Integers>,<n>) . . . . . factorization of an integer
##
IntegersOps.Factors := function ( Integers, n )
    return FactorsInt( n );
end;


#############################################################################
##
#F  FactorsInt( <n> ) . . . . . . . . . . . . . . prime factors of an integer
#F  FactorsRho(<n>,<inc>,<cluster>,<limit>) . . .  Pollards rho factorization
##
##  'FactorsInt' returns a list of prime factors of the integer <n>.
##
##  'FactorsInt' does trial divisions by the primes less than 1000 to  detect
##  all composites with a factor less than 1000 and primes less than 1000000.
##  After that it calls 'FactorsRho(<n>,1,16,8192)' to do the hard work.
##
##  'FactorsRho'  will  return a  list  of factors   and a list  of composite
##  number.   Usually  'FactorsInt'  factors  integers  with   prime  factors
##  $\<1000$ faster.     However  for   integers  with  no   factor  $\<1000$
##  'FactorsRho' will be faster.
##
##  'FactorsRho' uses Pollards $\rho$ method to factor the integer $n = p q$.
##  For a small simple example lets assume we want to factor $667 = 23 * 29$.
##  'FactorsRho' first calls 'IsPrimeInt' to avoid trying to factor a prime.
##
##  Then it uses the sequence defined by  $x_0=1, x_{i+1}=(x_i^2+1)$ mod $n$.
##  In our example this is $1, 2, 5, 26, 10, 101, 197, 124, 36, 630, .. $.
##
##  Modulo $p$ it takes on at most $p-1$ different values, thus it eventually
##  becomes recurrent, usually this happens after roughly $2 \sqrt{p}$ steps.
##  In our example modulo 23 we get $1, 2, 5, 3, 10, 9, 13, 9, 13, 9, .. $.
##
##  Thus there exist pairs $i, j$ such that $x_i = x_j$ mod $p$,  i.e.,  such
##  that $p$ divides $Gcd( n, x_j-x_i )$.  With a bit of luck no other factor
##  of $n$ divides $x_j - x_i$ so we find $p$ if we know such a pair.  In our
##  example $5, 7$ is the first pair, $x_7-x_5=23$, and $Gcd(667,23) = 23$.
##
##  Now it is too expensive to check all pairs, but there also must be  pairs
##  of the form $2^i-1, j$ with $3*2^{i-1} <= j < 4*2^{i-1}$.  In our example
##  $7, 13$ is the first such pair, $x_13-x_7=506$, and $Gcd(667,506) = 23$.
##
##  Thus by taking the gcds of $n$ and $x_j-x_i$ for such pairs, we will find
##  the factor $p$ after approximately $2 \sqrt{p} \<= 2 \sqrt^4{n}$ steps.
##
##  If $Gcd( n, x_j - x_i )$  is not a prime 'FactorsRho'  will  call  itself
##  recursivly with a different value for <inc>, i.e., it  will try to factor
##  the gcd using a different sequence $x_{i+1} = (x_i^2 + inc)$ mod $n$.
##
##  Since the gcd computations are by far the most time consuming part of the
##  algorithm  one can save time by  clustering differences and computing the
##  gcd  only every <cluster>  iteration.  This slightly increases the chance
##  that a gcd is composite, but reduces the runtime by a large amount.
##
##  Finally 'FactorsRho' accepts an argument <limit>  which is the number  of
##  iterations  performed by 'FactorsRho' before giving up. The default value
##  is  8192  which corresponds to a few minutes  while guaranteing that  all
##  prime factors less than $10^6$ and most less than $10^9$ are found.
##
##  Better descriptions of the algorithm and related topics can be found  in:
##  J. Pollard, A Monte Carlo Method for Factorization, BIT 15, 1975, 331-334
##  R. Brent, An Improved Monte Carlo Method for Fact., BIT 20, 1980, 176-184
##  D. Knuth, Seminumerical Algorithms  (TACP II),  AddiWesl,  1973,  369-371
##
FactorsRho := function ( n, inc, cluster, limit )
    local   i, sign,  factors,  composite,  x,  y,  k,  z,  g,  tmp;

    # make $n$ positive and handle trivial cases
    sign := 1;
    if n < 0  then sign := -sign;  n := -n;  fi;
    if n < 4  then return [ [ sign * n ], [] ];  fi;
    factors   := [];
    composite := [];
    while n mod 2 = 0  do Add( factors, 2 );  n := n / 2;  od;
    while n mod 3 = 0  do Add( factors, 3 );  n := n / 3;  od;
    if IsPrimeInt(n)  then Add( factors, n );  n := 1;  fi;

    # initialize $x_0$
    x := 1;  z := 1;  i := 0;

    # loop until we have factored $n$ completely or run out of patience
    while 1 < n  and 2^i <= limit  do

        # $y = x_{2^i-1}$
        y := x;  i := i + 1;

        # $x_{2^i}, .., x_{3*2^{i-1}-1}$ need not be compared to $x_{2^i-1}$
        for k  in [1..2^(i-1)]  do
            x := (x^2 + inc) mod n;
        od;

        # compare $x_{3*2^{i-1}}, .., x_{4*2^{i-1}-1}$ with $x_{2^i-1}$
        for k  in [1..2^(i-1)]  do
            x := (x^2 + inc) mod n;
            z := z * (x - y) mod n;

            # from time to time compute the gcd
            if k mod cluster = 0  then
                g := GcdInt( n, z );

                # if it is > 1 we have found a factor which need not be prime
                if g > 1  then
                    tmp := FactorsRho(g,inc+1,QuoInt(cluster+1,2),limit);
                    factors   := Concatenation( factors,   tmp[1] );
                    composite := Concatenation( composite, tmp[2] );
                              
                    n := n / g;
                    if IsPrimeInt(n)  then Add( factors, n );  n := 1;  fi;
                fi;
            fi;
        od;
    od;

    # add <n> to the list of composite numbers
    if 1 < n  then
        Add( composite, n );
    fi;

    # sort the list of factors and composite numbers and return it
    Sort(factors);
    Sort(composite);
    if 0 < Length(factors)  then
        factors[1] := sign * factors[1];
    else
        composite[1] := sign * composite[1];
    fi;
    return [ factors, composite ];

end;

FactorsInt := function ( n )
    local  sign,  factors,  p,  tmp;

    # make $n$ positive and handle trivial cases
    sign := 1;
    if n < 0  then sign := -sign;  n := -n;  fi;
    if n < 4  then return [ sign * n ];  fi;
    factors := [];

    # do trial divisions by the primes less than 1000
    # faster than anything fancier because $n$ mod <small int> is very fast
    for p  in Primes  do
        while n mod p = 0  do Add( factors, p );  n := n / p;  od;
        if n < (p+1)^2 and 1 < n  then Add(factors,n);  n := 1;  fi;
        if n = 1  then factors[1] := sign*factors[1];  return factors;  fi;
    od;

    # do trial divisions by known factors
    for p  in Primes2  do
        while n mod p = 0  do Add( factors, p );  n := n / p;  od;
        if n = 1  then factors[1] := sign*factors[1];  return factors;  fi;
    od;

    # handle perfect powers
    p := SmallestRootInt( n );
    if p < n  then
        while 1 < n  do
            Append( factors, FactorsInt(p) );
            n := n / p;
        od;
        Sort( factors );
        factors[1] := sign * factors[1];
        return factors;
    fi;

    # let 'FactorsRho' do the work
    tmp := FactorsRho( n, 1, 16, 8192 );
    if 0 < Length(tmp[2])  then
        Error( "sorry,  cannot factor ", tmp[2] );
    fi;
    factors := Concatenation( factors, tmp[1] );
    Sort( factors );
    factors[1] := sign * factors[1];
    return factors;
end;


#############################################################################
##
#F  DivisorsInt( <n> )  . . . . . . . . . . . . . . .  divisors of an integer
##
##  'DivisorsInt' returns a list of all divisors  of  the  integer  <n>.  The
##  list is sorted, so that it starts with 1 and  ends  with <n>.  We  define
##  that 'Divisors( -<n> ) = Divisors( <n> )'.
##
DivisorsSmall := [,[1],[1,2],[1,3],[1,2,4],[1,5],[1,2,3,6],[1,7]];
DivisorsInt := function ( n )
    local  divisors, factors, divs;

    # make <n> it nonnegative, handle trivial cases, and get prime factors
    if n < 0  then n := -n;  fi;
    if n = 0  then Error("DivisorsInt: <n> must not be 0");  fi;
    if n < 8  then return Copy( DivisorsSmall[n+1] );  fi;
    factors := FactorsInt( n );

    # recursive function to compute the divisors
    divs := function ( i, m )
        if Length(factors) < i     then return [ m ];
        elif m mod factors[i] = 0  then return divs(i+1,m*factors[i]);
        else return Concatenation( divs(i+1,m), divs(i+1,m*factors[i]) );
        fi;
    end;

    divisors := divs( 1, 1 );
    Sort( divisors );
    return divisors;
end;

DivisorsIntNonTriv := n -> let(d := DivisorsInt(n), d{[2..Length(d)-1]});

RPDivisorsIntNonTriv := n -> let(d := DivisorsInt(n), dd := d{[2..Length(d)-1]},
    Filtered(dd, x->Gcd(x, n/x)=1));

#############################################################################
##
#F DivisorPairs(x)
#F   returns a list of factorizations of x in two factors, 
#F   excluding [1,x] and [x,1]. Return type is
#F   a list of lists (each with 2 integers).
#F
DivisorPairs := function(x)
   local result;
   Constraint(IsInt(x));
   result := DivisorsInt(x);
   result := Sublist(result, [2..Length(result)-1]);
   result := List(result, i -> [i,x/i]);
   return result;
end;

#############################################################################
##
#F DivisorPairsRP(x)
#F   returns a list of factorizations of x in two relatively
#F   prime factors, excluding [1,x] and [x,1]. Return type is
#F   a list of lists (each with 2 integers).
#F
DivisorPairsRP := function(x)
   Constraint(IsInt(x));
   return Filtered( DivisorPairs(x), 
		    divpair -> Gcd(divpair[1], divpair[2]) = 1);
end;

#############################################################################
##
#F  Sigma( <n> )  . . . . . . . . . . . . . . . sum of divisors of an integer
##
##  'Sigma' returns the sum of the positive divisors of the integer <n>.
##
##  'Sigma' is a multiplicative arithmetic function, i.e., if $n$ and $m$ are
##  relative prime we have $\sigma(n m) = \sigma(n) \sigma(m)$.
##
Sigma := function ( n )
    local  sigma, p, q, k;

    # make <n> it nonnegative, handle trivial cases
    if n < 0  then n := -n;  fi;
    if n = 0  then Error("Sigma: <n> must not be 0");  fi;
    if n < 8  then return Sum(DivisorsSmall[n+1]);  fi;

    # loop over all prime $p$ factors of $n$
    sigma := 1;
    for p  in Set(FactorsInt(n))  do

        # compute $p^e$ and $k = 1+p+p^2+..p^e$
        q := p;  k := 1 + p;
        while n mod (q * p) = 0  do q := q * p;  k := k + q;  od;

        # combine with the value found so far
        sigma := sigma * k;
    od;

    return sigma;
end;


#############################################################################
##
#F  Tau( <n> )  . . . . . . . . . . . . . .  number of divisors of an integer
##
##  'Tau' returns the number of the positive divisors of the integer <n>.
##
##  'Tau' is a multiplicative arithmetic function, i.e., if $n$ and  $m$  are
##  relative prime we have $\tau(n m) = \tau(n) \tau(m)$.
##
Tau := function ( n )
    local  tau, p, q, k;

    # make <n> it nonnegative, handle trivial cases
    if n < 0  then n := -n;  fi;
    if n = 0  then Error("Tau: <n> must not be 0");  fi;
    if n < 8  then return Length(DivisorsSmall[n+1]);  fi;

    # loop over all prime factors $p$ of $n$
    tau := 1;
    for p  in Set(FactorsInt(n))  do

        # compute $p^e$ and $k = e+1$
        q := p;  k := 2;
        while n mod (q * p) = 0  do q := q * p;  k := k + 1;  od;

        # combine with the value found so far
        tau := tau * k;
    od;

    return tau;
end;


#############################################################################
##
#F  MoebiusMu( <n> )  . . . . . . . . . . . . . .  Moebius inversion function
##
##  'MoebiusMu'  computes the value  of  Moebius  inversion function for  the
##  integer <n>.   This  is 0 for  integers  which are not squarefree,  i.e.,
##  which are divided by a square $r^2$.  Otherwise it is 1 if <n> has a even
##  number and -1 if <n> has an odd number of prime factors.
##
MoebiusMu := function ( n )
    local  factors;

    if n < 0  then n := -n;  fi;
    if n = 0  then Error("MoebiusMu: <n> must be nonzero");  fi;
    if n = 1  then return 1;  fi;

    factors := FactorsInt( n );
    if factors <> Set( factors )  then return 0;  fi;
    return (-1) ^ Length(factors);
end;


#############################################################################
##
#F  CoefficientsQadic( <i>, <q> ) . . . . . .  <q>-adic representation of <i>
##
CoefficientsQadic := function( i, q )
    local   v;

    # represent the integer <i> as <q>-adic number
    v := [];
    while i > 0  do
        Add( v, RemInt( i, q ) );
        i := QuoInt( i, q );
    od;
    return v; 
end; 


#############################################################################
##
#F  IntegersOps.EuclideanDegree(<Integers>,<n>) . . . . . . . . absolut value
##
IntegersOps.EuclideanDegree := function ( Integers, n )
    if n < 0  then
        return -n;
    else
        return n;
    fi;
end;


#############################################################################
##
#F  IntegersOps.EuclideanRemainder(<Integers>,<n>,<m>)	. euclidean remainder
##
IntegersOps.EuclideanRemainder := function ( Integers, n, m )
    return RemInt( n, m );
end;


#############################################################################
##
#F  IntegersOps.EuclideanQuotient(<Integers>,<n>,<m>) . .  euclidean quotient
##
IntegersOps.EuclideanQuotient := function ( Integers, n, m )
    return QuoInt( n, m );
end;


#############################################################################
##
#F  IntegersOps.QuotientRemainder(<Integers>,<n>,<m>) . . . . . . quo and rem
##
IntegersOps.QuotientRemainder := function ( Integers, n, m )
    return [ QuoInt(n,m), RemInt(n,m) ];
end;


#############################################################################
##
#F  IntegersOps.QuotientMod(<Integers>,<r>,<s>,<m>)  quotient of two integers
#F                                                             modulo another
##
IntegersOps.QuotientMod := function ( Integers, r, s, m )
    if r mod GcdInt( s, m ) = 0  then
        return r/s mod m;
    else
        return false;
    fi;
end;


#############################################################################
##
#F  IntegersOps.PowerMod(<Integers>,<r>,<e>,<m>)  . . power of an integer mod
#F                                                                    another
#F  PowerModInt(<r>,<e>,<m>)  . . . . . . power of one integer modulo another
##
PowerModInt := function ( r, e, m )
    local   pow, f;

    # reduce r initially
    r := r mod m;

    # handle special case
    if e = 0  then
        return 1;
    fi;

    # if e is negative then invert n modulo m with Euclids algorithm
    if e < 0  then
        r := 1/r mod m;
        e := -e;
    fi;

    # now use the repeated squaring method (right-to-left)
    pow := 1;
    f := 2 ^ (LogInt( e, 2 ) + 1);
    while 1 < f  do
        pow := (pow * pow) mod m;
        f := QuoInt( f, 2 );
        if f <= e  then
            pow := (pow * r) mod m;
            e := e - f;
        fi;
    od;

    # return the power
    return pow;
end;

IntegersOps.PowerMod := function ( Integers, r, e, m )
    return PowerModInt( r, e, m );
end;


#############################################################################
##
#F  IntegersOps.Gcd(<Integers>,<n>,<m>) . . . . . . . . . gcd of two integers
#F  GcdInt(<n>,<m>) . . . . . . . . . . . . . . . . . . . gcd of two integers
##
IntegersOps.Gcd := function ( Integers, n, m )
    return GcdInt( n, m );
end;


#############################################################################
##
#F  IntegersOps.Lcm(<Integers>,<n>,<m>) . . least common multiple of integers
#F  LcmInt( <m>, <n> )  . . . . . . . . . . least common multiple of integers
##
IntegersOps.Lcm := function ( Integers, n, m )
    return LcmInt( n, m );
end;

LcmInt := function ( n, m )
    if m = 0  and n = 0  then
        return 0;
    else
        return AbsInt( m / GcdInt( m, n ) * n );
    fi;
end;


#############################################################################
##
#F  Gcdex( <m>, <n> ) . . . . . . . . . . greatest common divisor of integers
##
Gcdex := function ( m, n )
    local   f, g, h, fm, gm, hm, q;
    if 0 <= m  then f:=m; fm:=1; else f:=-m; fm:=-1; fi;
    if 0 <= n  then g:=n; gm:=0; else g:=-n; gm:=0;  fi;
    while g <> 0  do
        q := QuoInt( f, g );
        h := g;          hm := gm;
        g := f - q * g;  gm := fm - q * gm;
        f := h;          fm := hm;
    od;
    if n = 0  then
        return rec( gcd := f, coeff1 := fm, coeff2 := 0,
                              coeff3 := gm, coeff4 := 1 );
    else
        return rec( gcd := f, coeff1 := fm, coeff2 := (f - fm * m) / n,
                              coeff3 := gm, coeff4 := (0 - gm * m) / n );
    fi;
end;


#############################################################################
##
#F  Int( <obj> )  . . . . . . . . . . . . . . . . . . . convert to an integer
##
Int := function ( obj )
    if IsInt( obj )  then
        return obj;
    elif IsRat( obj )  then
        return QuoInt( Numerator( obj ), Denominator( obj ) );
    elif IsFFE( obj )  then
        return IntFFE( obj );
    else
        Error("<obj> must be rational or a finite field element");
    fi;
end;


#############################################################################
##
#F  AbsInt( <n> ) . . . . . . . . . . . . . . .  absolute value of an integer
##
AbsInt := function ( n )
    if 0 <= n  then return  n;
    else            return -n;
    fi;
end;


#############################################################################
##
#F  SignInt( <n> )  . . . . . . . . . . . . . . . . . . .  sign of an integer
##
SignInt := function ( n )
    if   0 =  n  then
        return 0;
    elif 0 <= n  then
        return 1;
    else
        return -1;
    fi;
end;


#############################################################################
##
#F  ChineseRem( <moduli>, <residues> )  . . . . . . . . . . chinese remainder
##
ChineseRem := function ( moduli, residues )
    local   i, c, l, g;

    # combine the residues modulo the moduli
    i := 1;
    c := residues[1];
    l := moduli[1];
    while i < Length(moduli)  do
        i := i + 1;
        g := Gcdex( l, moduli[i] );
        if g.gcd <> 1  and (residues[i]-c) mod g.gcd <> 0  then
            Error("the residues must be equal modulo ",g.gcd);
        fi;
        c := l * (((residues[i]-c) / g.gcd * g.coeff1) mod moduli[i]) + c;
        l := moduli[i] / g.gcd * l;
    od;

    # reduce c into the range [0..l-1]
    c := c mod l;
    return c;
end;


#############################################################################
##
#F  LogInt( <n>, <base> ) . . . . . . . . . . . . . . logarithm of an integer
##
LogInt := function ( n, base )
    local   log;

    # check arguments
    if n    <= 0  then Error("<n> must be positive");  fi;
    if base <= 1  then Error("<base> must be greater than 1");  fi;

    # 'log(b)' returns $log_b(n)$ and divides 'n' by 'b^log(b)'
    log := function ( b )
        local   i;
        if b > n  then return 0;  fi;
        i := log( b^2 );
        if b > n  then return 2 * i;
        else  n := QuoInt( n, b );  return 2 * i + 1;  fi;
    end;

    return log( base );
end;

#############################################################################
##
#F  CeilLog2( <n> ) . . . . . . . . . .base 2 logarithm rounded up to integer
##
CeilLog2  := x -> IntDouble(d_ceil(d_log(x)/d_log(2)));
FloorLog2 := x -> IntDouble(d_floor(d_log(x)/d_log(2)));

#############################################################################
##
#F  RootInt( <n>, <k> ) . . . . . . . . . . . . . . . . .  root of an integer
##
RootInt := function ( arg )
    local   n, k, r, s, t;

    # get the arguments
    if   Length(arg) = 1  then n := arg[1];  k := 2;
    elif Length(arg) = 2  then n := arg[1];  k := arg[2];
    else Error("usage: 'Root( <n> )' or 'Root( <n>, <k> )'");
    fi;

    # check the arguments and handle trivial cases
    if  k <= 0                  then Error("<k> must be positive");
    elif k = 1                  then return n;
    elif n < 0 and k mod 2 = 0  then Error("<n> must be positive");
    elif n < 0 and k mod 2 = 1  then return -RootInt( -n, k );
    elif n = 0                  then return 0;
    elif n <= k                 then return 1;
    fi;

    # r is the first approximation, s the second, we need: root <= s < r
    r := n;  s := 2^( QuoInt( LogInt(n,2), k ) + 1 ) - 1;

    # do Newton iterations until the approximations stop decreasing
    while s < r  do
        r := s;  t := r^(k-1);  s := QuoInt( n + (k-1)*r*t, k*t );
    od;

    # and thats the integer part of the root
    return r;
end;


#############################################################################
##
#F  SmallestRootInt( <n> )  . . . . . . . . . . . smallest root of an integer
##
SmallestRootInt := function ( n )
    local   k, r, s, p, l, q;

    # check the argument
    if   n > 0  then k := 2;  s :=  1;
    elif n < 0  then k := 3;  s := -1;  n := -n;
    else return 0;
    fi;

    # exclude small divisors, and thereby large exponents
    if n mod 2 = 0  then
        p := 2;
    else
        p := 3;  while p < 100  and n mod p <> 0  do p := p+2;  od;
    fi;
    l := LogInt( n, p );

    # loop over the possible prime divisors of exponents
    # use Euler's criterion to cast out impossible ones
    while k <= l  do
        q := 2*k+1;  while not IsPrimeInt(q)  do q := q+2*k;  od;
        if PowerModInt( n, (q-1)/k, q ) <= 1  then
            r := RootInt( n, k );
            if r ^ k = n  then
                n := r;
                l := QuoInt( l, k );
            else
                k := NextPrimeInt( k );
            fi;
        else
            k := NextPrimeInt( k );
        fi;
    od;

    return s * n;
end;


#############################################################################
##
#F  IntegersOps.AsGroup(<Integers>) . . . . . . . . . .  view the integers as
#F                                                       multiplicative group
##
IntegersOps.AsGroup := function ( Integers )
    Error("sorry, Z is not finitely generated as multiplicative group");
end;


#############################################################################
##
#F  IntegersOps.AsAdditiveGroup(<Integers>) . . . . . .  view the integers as
#F                                                             additive group
##
#N  14-Oct-91 martin this should be
#N  IntegersAsAddtiveGroupOps := Copy( AdditveGroupOps );
##
IntegersAsAdditiveGroupOps := Copy( DomainOps );
IntegersAsAdditiveGroupOps.\in   := IntegersOps.\in;
IntegersAsAdditiveGroupOps.Random := IntegersOps.Random;

IntegersOps.AsAdditiveGroup := function ( Integers )

    return rec(

        isDomain                := true,
        isAdditiveGroup         := true,

        generators              := [ 1 ],
        zero                    := 0,

        isFinite                := false,
        size                    := "infinity",
        isCyclic                := true,

        operations              := IntegersAsAdditiveGroupOps
    );

end;


IsPosInt := x->IsInt(x) and x > 0;
IsNegInt := x->IsInt(x) and x < 0;
IsPosInt0 := x->IsInt(x) and x >= 0;
IsNegInt0 := x->IsInt(x) and x <= 0;

Is2Power := n -> Cond(IsInt(n), 2^Log2Int(n)=n, 
                      spiral.code.IsValue(n), 2^Log2Int(n.v)=n.v,
		      false);
    
IsIntPower := (n, exp) -> exp^LogInt(n, exp) = n;

